Real-Time Dynamics in Quantum Impurity Models with Diagrammatic Monte Carlo 
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We extend the recently developed real-time Diagrammatic Monte Carlo method, in its hybridiza- 
tion expansion formulation, to the full Kadanoff-Baym-Keldysh contour. This allows us to study 
real-time dynamics in correlated impurity models starting from an arbitrary, even interacting, initial 
density matrix. As a proof of concept we apply the algorithm to study the non equilibrium dynamics 
after a local quantum quench in the Anderson Impurity Model. Being a completely general approach 
to real-time dynamics in quantum impurity models it can be used as a solver for Non Equilibrium 
Dynamical Mean Field Theory. 
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I. INTRODUCTION 



The understanding of real-time dynamics in strongly 
correlated quantum systems represents a major challenge 
in modern condensed matter physics, due to the rapid 
experimental advances in probing physical responses di- 
rectly in the time domain^. Time resolved pump-probe 
experiments with femtosecond resolution, for example, 
have been recently moving from the realm of atoms 
and molecules 2 - to that of strongly correlated bulk ma- 
terials, such as Mott Insulators 3 * and High- Temperature 
Superconductors 4 -, thus opening the way toward a full 
characterization of their non equilibrium properties. Sim- 
ilarly, the emerging field of cold atomic quantum gases 
represents a natural laboratory where the dynamics of 
almost isolated quantum systems can be probed in real- 
time. Several experiments have been performed in this 
direction to measure, for instance, the dynamics after a 
quantum quench in bosonic gases 5 , the lifetime of dou- 
blons in strongly interacting Fermi systems^ or the onset 
of superexchange interactions between ultra cold atoms 
in optical lattices^. Finally, we note that even at the 
nanoscale level the experimental frontiers are rapidly 
moving in the direction of time-resolved techiniques to 
detect charge transport by counting individual electrons 
while tunneling across correlated nanostructures such as 
semiconducting quantum dots^r— . 

From the theoretical point of view, the field of real- 
time dynamics in strongly correlated systems is still at 
its infancy. In this perspective Quantum Impurity (QI) 
models represent the ideal playground to test and de- 
velop new methods. These models consist of a small 
quantum system with few interacting degrees of freedom, 
the impurity, coupled via hybridization to a reservoir of 
fermionic excitations. QI models therefore represent, by 
construction, the natural framework to study quantum 
transport through nanocontacts. While the equilibrium 
physics of these nutshell strongly correlated systems can 
be studied with a wide range of powerful numerical and 
analitical tools, their non equilibrium real-time dynamics 
is still challenging. The reason for this gap is mainly due 



to the fact that most of the theoretical tools which has 
been developed in the last thirty years to solve quantum 
impurity models in equilibrium, most notably Numerical 
Renormalization Group 11 ' 12 (NRG), can not be directly 
applied to the out of equilibrium case. This has trig- 
gered a large amount of theoretical works among which 
we mention the time dependent extensions of NRGi 3 - and 
DMRG— li 5 -, the ISPI method^ and the Flow equation 
approach^. 

Beside their relevance for nanoscience, QI models have 
also been emerging in the last two decades as the 
paradigm to understand strong correlation phenomena in 
bulk lattice models within the so called Dynamical Mean 
Field Theory^. The extension of this powerful non per- 
turbative technique to the non equilibrium domai n 19 : 20 
makes the development of a generic numerically exact ap- 
proach to real-time dynamics in quantum impurities an 
even more urgent issue. 

Recently, a new generation of Diagrammatic Monte 
Carlo (diagMC) algorithms have been developed to solve 
equilibrium quantum impurity models. These are based 
on a stochastic sampling of the partition function written 
as an imaginary-time diagrammatic expansion around 
weak or strong coupling values of the interaction acting 
on the impurit y 21 ' 22 . An extension to the out of equilib- 
rium case, namely to the real-time domain, has been pro- 
posed only very recently. The general idea of Real-Time 
diagMC methods is to start from a given initial density 
matrix, describing the fermionic reservoir and the impu- 
rity, and to compute the dynamics of any observable of 
the system by sampling the real-time evolution operator 
written as a diagrammatic expansion along the Keldysh 
contour. Both the weak and the strong coupling expan- 
sion have been proposed and applied to the local polaron 
proble m 23 : 24 and to the Anderson impurity model 25 . By 
construction, these approaches rely on a rather specific 
assumption on the initial density matrix, which has to 
describe either a non-interacting impurity model or an 
impurity decoupled from the reservoir. Going beyond 
this assumption, which is technical rather than funda- 
mental, is the main motivation for this work. To this 
extent we merge together the imaginary time and the 
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real-time methods thus developing a completely general 
Diagrammatic Monte Carlo algorithm on the full three- 
branches Kadanoff-Baym-Keldysh contou r 26 ! 27 . This al- 
lows us to deal with the most generic non equilibrium 
set-up, namely an interacting quantum impurity model 
described at time t = by some thermal density matrix, 
driven out of equilibrium for time t > by a generic, 
possibly time dependent, perturbation. As a by product, 
we obtain a general real-time solver for quantum impu- 
rity models which can be therefore used to solve Non 
Equilibrium Dynamical Mean Field Theory. Even in this 
context, the possibility of studying the non equilibrium 
real-time dynamics starting from interacting initial states 
look particularly intriguing in the light of first DMFT 
results on quantum quenches in the Hubbard model 20 
which focus so far only on quenches starting from a non 
interacting initial state. 

As a first application of this algorithm we study the 
non equilibrium real-time dynamics in the Anderson Im- 
purity model after a local quantum quench. 

The paper is organized as follows. In section II we 
introduce the formalism. In section III we derive the 
hybridization expansion on the Kadanoff-Baym-Keldysh 
contour while section IV is devoted to the description of 
the diagMC algorithm. Section V describes the applica- 
tion to the Anderson Impurity Model while section VI is 
for conclusions and perspectives. 

II. NON EQUILIBRIUM DYNAMICS IN 
QUANTUM IMPURITY MODELS 

The aim of this section is to set-up the proper frame- 
work to study non equilibrium real-time dynamics in 
quantum impurity models. To this purpose, we consider 
a set of discrete electronic levels, the "impurity", with 
creation operator c\ labeled by an integer a — 1, . . .J\f 
which may include both spin and orbital degrees of free- 
dom. These levels are coupled to one or more baths of 
free fermions with momentum k and creation operator 
The generic hamiltonian of a QI model reads 

= £ k fLha +nj oc [4,C ] 

k a 

+ J2( V ^fL^a + h.C.) , (1) 
k a 

where the first term describes the continuum of fermionic 
excitations, the local hamiltonian T~L[oc \_ c \n c a\ generally 
contains many-body interactions for electrons on the im- 
purity, while the last term is the hybridization which cou- 
ples the impurity and the bath and it is assumed here, 
for the sake of simplicity, to be diagonal in the a index. 

Since we are interested in studying non equilibrium 
dynamics of model (fl}, we have to specify an initial 
condition as well as a protocol to drive this system out of 
equilibrium. Following general ideas of non equilibrium 
many body theor y 26 ! 27 , we imagine to prepare the system 



at t = in a thermal state of namely we choose the 
Boltzmann distribution as initial density matrix 

e -pH- 

p(t = 0) = Peq = ^^, Z = Tve-' 3H -, (2) 

and then, for t > 0, let the system evolve under the action 
of a new hamiltonian 

H{t) =«_ + V(f) , t>0 (3) 

Choosing the initial density matrix as the thermal one 
gives access to the response of a correlated quantum im- 
purity model to external fields. As we shall discuss in 
what follows, Eq. @ represents the main point where our 
approach differs from previous implementations of the 
real-time diagrammatic Monte Carlo method 2 ^—. For 
what concerns the driving protocol, namely the nature 
of the external perturbation, there are actually several 
ways to push a quantum impurity model out of equilib- 
rium. In this work we shall focus on the simplest one, 
namely a quantum quench experiment, but the method 
allows to address even more general time dependent out 
of equilibrium problems. 

In a quantum quench, one imagines to prepare the sys- 
tem, at t — 0, in a given state of some initial hamiltonian 
(H- in the case of our interest) and then, for t > 0, to 
suddenly change some of its parameters letting evolve the 
system under the unitary action of a new Hamiltonian 
H+. Such a protocol therefore represents the simplest 
example of time-dependent problem where the variation 
in time is step-like 

H{t) = H- + 9{t)8U, 5H = n + -H-. (4) 

The sudden quench injects energy into the system and 
leads to a relaxation dynamics towards a new steady 
state, provided the perturbation STi is not a conserved 
quantity of The main task is therefore to com- 

pute quantum averages with the full density matrix p(t) 
evolved in real-time 

(0(t))=Tr[p(t)0\ =Tr[ Peq uHt)OU(t)] . (5) 

where the trace has to be taken over the bath and the 
impurity degrees of freedom, while U (t) and W (t) are, 
respectively, the unitary operator generating the dynam- 
ics and its hermitian conjugate. In the specific case of a 
time independent hamiltonian, as we have for t > see 
Eq. (H|), these operators read 

U(t) = e- in + t U^(t) = e in + t . (6) 

To proceed further, it is convenient to specify the nature 
of the perturbation 5 % induced by the quantum quench. 
To keep the discussion as general as possible we write the 
hamiltonian of the system for t > as 

k a 

+ JL{ V ^afL^ + h.c) , (7) 
k a 
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namely we allow for an abrupt change of all the pa- 
rameters entering in the hamiltonian ([1]), in such a way 
that different kind of non equilibrium phenomena can be 
treated within the present approach. One can, for ex- 
ample, study the dynamics after a local quantum quench 
acting only on the impurity degrees of freedom, or con- 
sidering a global change in the hamiltonian, acting on 
the hybridization term or even on the conduction band. 
These global quantum quenches are particularly relevant 
for studying dc transport in quantum dots and will be 
the subject of a forthcoming publication. 

Once we have specified the structure of the hamil- 
tonian after the quench, we can perfom the hybridiza- 
tion expansion in complete analogy to what has been 
done previously in the case of a pure Keldysh real-time 
algorithm 23-25 with, nonetheless, an important differ- 
ence. Indeed, as it appears clearly from the above formu- 
lation (see Eq.©, not only the real-time evolution but 
also the "thermal" one is governed by the full hamiltonian 
7i± , involving both the impurity and the bath degrees of 
freedom coupled one to each other. This allows us to 
overcome the limitation of the pure Keldysh approach 
which relies on the assumption of an initially decoupled 
density matrix and gives us the possibility to treat ar- 
bitrary initial conditions. In the next sections we will 
describe the Diagrammatic Monte Carlo algorithm used 
to compute the real-time average in Eq. ([5]). 



III. HYBRIDIZATION EXPANSION ON THE 
KADANOFF-BAYM-KELDYSH CONTOUR 



In order to study the non equilibrium real-time dynam- 
ics of quantum impurity models starting from a generic 
initial density matrix, we formulate the diagrammatic 
monte carlo algorithm (diagMC), in its hybridization ex- 
pansion version, on the Kadanoff-Baym-Keldysh contour 
made by the usual imaginary time axis and the real-time 
Keldysh contour. As we are going to show, this struc- 
ture naturally emerges from the definition of real-time 
quantum averages given in equation ([S]) . To proceed fur- 
ther, we introduce a dynamical time-dependent partition 
function for the QI model which is defined as 



the real-time operator and its hermitian conjugate we get 



Z(t,p) = Tr [e-^ Ua U^ {t)U{t)\ 



(8) 



We note that this quantity does not actually depend on 
time t since, by construction, the evolution is unitary, 
nevertheless it represents the natural quantity to derive 
the hybridization expansion. As it will appear more 
clearly later on, Z (t, /3) can be seen a dynamical gen- 
erating functional of the Monte Carlo weights needed to 
compute any quantum average in real-time. The basis of 
any continuous-time diagMC algorithm is the expression 
of evolution operators as time-ordered exponentials. For 



U (t) = T exp 



Eft (t) = T exp (ij dtU 



i I dtU+{t) 
o 

/ 

(*) 



(9) 
(10) 



where T (T) is the time ordering (anti-time ordering) 
operator, whose action is order a string of real-time 
fermionic operators according to their time arguments, 
placing to the left the operators at later (earlier) times, 
with an overall plus or minus sign according, respectively, 
to the parity of the number of fermionic exchanges needed 
to move the string from the original to the final ordering. 
Using the well known properties of the equilibrium den- 
sity matrix ([2]) we can write also the Boltzmann weight 
as an imaginary time evolution along the path [— ij3, 0] 



-pv.- 



T exp £ dT K- (-ir)j 



(11) 



= T^ exp —i 



-i/3 



dtU- (t) =U(-if3) 



where T-| is an imaginary-time-ordering operator defined 
in complete analogy with T. 

Inserting these expressions in the dynamical partition 
function previously introduced, we get 



Z(t) 



= Tr 
= Tr 



T n e 



if° ifj dtu-(t) ^ el / o ' dtn+{t) Te if° dtu+(t) 



(12) 



where, in the second line, C is the Kadanoff-Baym- 
Keldysh contour plotted in figure [TJ which is made of 
three branches B,, i — 1,2,3, being respectively the up- 
per real-time branch, the lower one and the imaginary- 
time branch. Hereafter the time argument t in Eq. (fT2|) 
is assumed to live on such a contour, unless differently 
specified. Time ordering along C, enforced by operator 
Tc , acts similarly to the standard Keldysh time-ordering, 
placing operators with later time on the left, namely 



T C (A (t 1 )A{t 2 )) = 



A(ti)A(i 2 ) h>t 2 



-A(h)A(h) h < t 2 



(13) 



G ,C\ 



where > (<) means greater (lesser) on the contour C. By 
the definition of the partition function Z(t), it follows 
that 



Tr 



T] ®T0T. 



The integral along the contour is defined as 



dt = 



c 



dt - 



i/3 



dt 



dt, 



(14) 



(15) 
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Figure 1: Kadanoff-Baym-Keldysh contour which starts at 
time on the first branch, runs up to time t — t* then back 
from i, to t = and finally along the imaginary axis from 
i = to t = -ift. 



while, along the contour, the Hamiltonian entering in 
Eq. (T5J is 



H(t) 



teB 1 ,B 2 
t£B 3 



H- 



(16) 



Once we have written the partition function Z (t, ft) as a 
time ordered exponential, we can treat different terms in 
the contour hamiltonian as c-numbers. Therefore we can 
write Eq. (fl~2"j) as 



Z(t,ft)=Tr 



if dtn (t)+H hyb (t) 



(17) 



where we have explicitly indicated the hybridization 
Hamiltonian Whyb (t) , defined as 



n hyb (t) = ^2(v ka (t) fljt) Ca (t) 

k a 



h.c 



with a time dependent hybridization Vk a (t) — V^ a de- 
pending on the position of time t along the contour C. 
To proceed further, it is convenient to introduce the bath 
operators at the impurity site, defined as 



ft(t) = Y,v ka (t)fl a (t) 



k 



(19) 



which enter the hybridization Hamiltonian Eq. (|18l) . 
Then we formally expand Z (t, ft) in power of the hy- 
bridization Hamiltonian (fT8|) and trace out, at any order 
in the expansion, separately the bath and the impurity 
degrees of freedom which are completely decoupled in the 
absence of Tihyb- Let us define n a & {n a b) as the number 
of creation(annihilation) impurity operators, in the fol- 
lowing also called kinks, with flavour a and in branch b. 
These integers run in general between zero and infinity, 
the only constraint is that the total number of creation 
and annihilation kinks for each flavour a, say k a and k a , 
have to be equal due to total particles conservation. This 
introduces the following constraint 



k a nab = k a = ^2 n a b a = l,...,M 



(18) The resulting expansion for the partition function reads 



z{t,ft) = nn <n ab ,n ab ) n /< /^n( T c/(*?)/ t (*?)---/(*L)/ t (*i)>^ 

1 s - _ - i J J a 

x ( T C ft (ct (tl) c (t?) . . . c t (tgj c (t*J) ) local , (20) 



a—l b—1 n a b.n a 



where the contour integrals must be constrained to the 
time ordered regions 

c c c - c _ c c ~ , . 

*?><§>...> 3. , ti>n> ...>t% a , (21) 

while s(n a b, n a b) includes all the signs(phases) coming 
from the real-time evolution operators as well as from 
the integration along the imaginary time axis 

S(n ab , n ab ) = (_*)«■»+»•» i fia2+na 2 (.j)*..^.. . ( 22 ) 

By using our definition of the total number of kinks k we 
can write these factors as 

s(n ab , n ab ) = (-!)*» (-1)"-+"- (-ifas+n a3 (23) 



Concerning the trace over the bath degrees of freedom 
with flavour a, this can be written by using Wick's the- 
orem as the determinant of a k a X k a matrix A a 

( T C f (tj) P (t$) . . . / (tlj /t (^J) tath - detA a (24) 

whose entries A"- are the contour-ordered hybridization 
functions defined as 

= iA c (t«,t«) = 

= £ ^ a (*?) ( T C / k (t?) /t ) F k a (t?) (25) 
k 

the average being taken over the bath degrees of free- 
dom. This function is defined along the Kadanoff-Baym- 
Keldysh contour and therefore naturally acquires the 
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structure of a 3 x 3 matrix in the branch spaced, as 
we discuss in the appendix |XJ 

Unlike what happens for the bath, the trace over the 
local hilbert space cannot be written in terms of sin- 
gle particle contractions since in general Wick's theorem 
does not hold for an interacting impurity. This consid- 
eration holds regardless of the specific form of the local 
hamiltonian H.i oc and it is also true in equilibrium. It is 
therefore clear that the evaluation of the local term in 
Eq. (f2T)|) for a given configuration of kinks is the compu- 
tational bottlneck of the algorithm, expecially in the case 
where multi-orbital interactions are considered. We note 
that such a kind of expressions arise also in NCA-kind 
of approaches to QI models 29 . To efficiently evaluate 
this local term we follo w 30 ' 31 and take advantage of the 
reduced hilbert space of the impurity to write the multi- 
point correlation function of Eq. (I20[) in the basis of the 
local eigenstates. The trace then reduces to multiplying 
matrices sandwiched by local evolution operators. If we 



define a global time ordering along the contour such that 

t\ > *2 > • • • > hN , 

where N — J2 a — J2 a ka> t nen the local trace can be 
written as 

(T C Ua (*?) C (t?) ■ • • Ct (tgj C (tgj) ) local = 

s To Tr [pioc X 1 (h) X 2 (h) ...X N (t N ) ] (26) 

where we have introduced an extra sign st c due to time 
ordering. In the previous equation the X's are creation or 
annihilation operators (depending on the time ordering) 
evolved in time with the local hamiltonian 

Xi (ti) = e iHloatl Xi e - in ^ti l = l,...,2N, (27) 

Combining all the above results, we write the hybridiza- 
tion expansion for the dynamical partition function as 



Z(t,j3) 



AT 3 

nn 



s(n a 



b i ^ab I 



H / dt i / <II Det t A<I ] ToftlPiocXiih) ...X 2N (t 2N )} 



(28) 



-1 b— 1 n a b,n a 



It is worth noticing that this expression represents an ex- 
act result for the partition function of the original quan- 
tum impurity model. As we are going to show in the 
next section, the goal of the Diagrammatic Monte Carlo 
method is to sum up stochastically the hybridization ex- 
pansion using a Metropolis algorithm. 

A. Effective Action Formulation and 
Non-Equilibrium DMFT 

An important outcome of the previous calculation has 
been to obtain an exact and closed expression for the dy- 
namical partition function Z (t, /3) of the quantum impu- 
rity model, written as a functional of the contour-ordered 
hybridization Ac(t,t'). This result holds in general. 
Indeed, right the same expression for Z (t, (3) is recov- 
ered within a path integral formulation of the problem. 
Following Kamene\i22, we define the dynamical partition 
function Z (t) as a path integral over the fermionic co- 
herent states c (t) , c (t) defined along the three branch 
contour 

Z(t)= Vc a Vc a e iS 'ff^-^ (29) 

a . 

The effective action describing the real-time dynamics of 
the impurity model can be written generally as 

iS ef f=i Sloe + f dt [ dt'J^ C " (0 iA C (*> 5 a (*) 

J c J c a 

(30) 



where Si oc is the local-in-time impurity action while the 
quadratic part is defined in terms of the contour-ordered 
hybridization function iAc(t,t') which takes into ac- 
count the coupling to the bath. It is worth noticing here 
that, by construction, time arguments lie along the three- 
branch contour C, while the integrals are defined as in Eq. 
(fT"5)) . As a consequence, the contour-ordered hybridiza- 
tion function naturally acquires a 3 x 3 matrix structure 
in the Keldysh-Matsubara space, 



iA c (*,*') -> iA aP (t,t') a, = 1,2,3 (31) 



where the 2x2 block with a, (3 ^ 3 is the Keldysh 
subspace, the last diagonal element a — (3 — 3 is the 
Matsubara sector while the remaining off-diagonal terms 
are mixed hybridization functions describing the effect 
of the initial condition. We want to show now how is 
possible within this framework to recover the result for 
hybridization expansion we previously derived. The idea 
is to formally expand the effective action in power of the 
contour-ordered hybridization iAc and use the defini- 
tion of path-integral as the contour time-ordered average 
of field-operatorsS 2 . to get for the partition function Z (t) 
the following expansion 
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a fe a a ' ^ C ^ C a 

, (Ml 



where the average is taken over the initial local den- 
sity matrix. To proceed, we symmetrize the integrand 
with respect to the k a \ permutations of creation times 
(t", <2, • • • } resulting into an extra l/k a \ and a deter- 
minantal combination of the k a hybridization functions, 
where the correct signs to build the determinant are pro- 
vided by the contour time ordered local trace. The fc Q -th 
order term in the expansion therefore reads 



n, 



(-ir° 

(TcU 



Ic dt t 



dt ka Det [A a ] 



(*!))• (33) 



Now the integrand is fully symmetric under permutations 
and we can take advantage of this fact to reduce the size 
of the integration domain, namely we choose among the 
k a \ possible contour time-orderings for the creation and 
annihilation times with flavour a the following ones 



c c 
t a 1 >t%> 



_ c _ c 

t? > ts > 



C ~ 



(34) 



As anticipated, the final result we get for the dynamical 
partion function Z (t) coincides with the one quoted in 
Eq. ([28]). 

These considerations are relevant for studying quan- 
tum quenches and real-time dynamics of correlated 
lattice models within Dynamical Mean Field Theory 
(DMFT). In this case, as one can show explicitly, the 
full non equilibrium many body problem is mapped, in 
the limit of infinite lattice coordination, onto a quantum 
impurity model coupled to a non- equilibrium bath and 
subject to a self-consistency condition. The dynamical 
partion function of this effective non equilibrium problem 
acquires exactly the same form as in Eq. (|29|) , with an un- 
known contour ordered hybridization function iAc (t, t'), 
which generally lacks time traslational invariance. This 
fact makes uneffective most of the conventional impu- 
rity solvers used in equilibrium DMFT, which rely on 
a time independent hamiltonian formulation of the effec- 
tive problem, thus suggesting diagrammatic Monte Carlo 
method as a natural candidate to solve Non-Equilibrium 
Dynamical Mean Field equations. 



which is obtained by a proper extension of the graphical 
representation introduced by Werner and coworker o 22 i 30 
in the context of the imaginary-time diagMC algorithm. 

As it can be immediately read out from Eq. (1251) . 
the dynamical partition function can be written as a 
weighted sum over configurations C made by diagrams 
on the Kadanoff-Baym-Keldysh contour 



(35) 



A given configuration contains, for each channel a — 
1, ... N, a total of 2k a vertices occurring at times {tf , tf} 
on the contour, with i = 1, . . . k a . Half of these vertices 
represent an impurity creation operator c\ (tf) while the 
other half stems for an impurity annihilation operator 
c a (tf), both of them being evolved in time with the local 
hamiltonian "Hi oc . All together we have 2j2 a impu- 
rity operators which we store in such a way to always 
preserve global time ordering along the contour. In sum- 
mary a typical configuration reads 



a = 1,2,--- ,JV 

k a = 0, . . . 00 

C — I c c c 

*1 > *2 > • • • > *fe a 

. c ~ c c ~ a 
*?>*£>...>« 



(36) 



An example of such a configuration is shown in figure [5] 
The weight W (C) associated to each configuration C can 

O • „ 



t -o- 
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Figure 2: An example of configuration C sampled by diagMC 
algorithm. 



IV. DIAGRAMMATIC MONTE CARLO 

Diagrammatic Monte Carlo is a numerical algorithm 
for sampling infinite series of multiple integrals, such as 
those arising in any perturbative expansion 33 . As it is 
well known in many body theory, these expansions often 
admit a diagrammatic representation 3 ^, even in out-of- 
equilibrium situations. This is true also for the hybridiza- 
tion expansion of section lnTl as we are now going to show, 



be read directly from the hybridization expansion, see 
Eq. (|28p. For later convenience we define it as 

W (C) = Det [C] sign [C] Tr ioc [C] , (37) 

where sign [C] includes all the signs (phases) coming from 
the evolution as well as from the time ordering, while the 
trace over the configuration reads 

Tr /oc [C] = Tr [p loc X l (t t ) X 2 (i 2 ) . . . X 2N [t 2N ) ] , (38) 
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so that by definition we get for Z {€) 



(39) 



The same weights W (C) are also required for evaluating 
the average of any local operator O acting on the im- 
purity Hilber space. Indeed, if consider the definition of 
quantum averages given in Eq. ([S]), namely 



(0(t)}=Tr [p eq U^ (t) OU{t)} 



(40) 



and perform the hybridization expansion of the previous 
section, we find that 



(o(t)) = 



C{C)W(C) 



(41) 



where the estimator of local operator has been defined as 



0(C) 



Tr[p; oe Xi(ti) ...O... X 2N (t 2N ) 

Tr [pi oc X\ . . . X 2N (*2Ar)] 



(42) 



Once the real-time average of a local operator is writ- 
ten as in Eq. (|41| , it would be natural natural to sample 
it using a Monte Carlo method, namely generating a ran- 
dom walk in the configuration space which visit configu- 
rations C with probability P (C) = W (C) / J2c W 

When trying to implement this idea in the context of 
real-time quantum dynamics the problem one have to 
face is that the weight W (C) is in general a complex 
number. In the specific case of interest this is not only 
due to the explicit "i-factors" coming from the real-time 
evoultion and entering sign[C] but also to the fact that 
the contour ordered bath, defined in Eq (|23|) and enter- 
ing the determinants, is indeed a complex function of its 
time arguments (see Appendix [A|. The simplest way to 
circumvent this problem is to sample the absolute value 
of the weight, \W (C) |, while including the phase r\ (C), 
defined as 



W(C) 
\W(C)\ 



(43) 



in the Monte Carlo estimator. While this approach al- 
lows for a straightforward implementation, it becomes 
problematic when the average phase goes to zero. In 
this respect, we note that more refined but computation- 
ally expensive techniques, based on sampling blocks of 
configurations at fixed sign, has been developed in recent 
years to cope with this dynamical sign problem 3 ^ . There- 
fore a possible future direction could be to merge them 
with present implementations of diagMC method to see 
if a compromise between efficiency and accurancy can be 
found. For the time being we avoid this route, sampling 
directly the absolute value of the weight. 



A. Metropolis Algorithm 

The standard technique to generate configurations 
with a given probability (in the case of our interest 



Poo (C) = \W (C) \/J2c> \ W i C ') I) is t0 build U P a Markov 
chain, namely a stochastic process which describes the 
evolution of the probability to visit configuration C after 
n steps. 



P(C,n) = Proba(C(n) = C) 



(44) 



This Markov chain is fully determined once we assign the 
conditional probability S [C \ C] to be in C at step n + 1 
being in C at step n. Indeed, this is the quantity entering 
the master equation 



P(C',n+l) = S[C'\C] P(C,n) 



(45) 



Sufficient conditions for this master equation to reach, af- 
ter waiting a proper equilibration time, the desired prob- 
ability Pqo (C) is that the matrix S [C'\C] is ergodic and 
satisfies the so called detailed balance condition. Ergod- 
icity means that it has to be possible to reach any con- 
figuration C from any other configuration C in a finite 
number of steps, while detailed balance means that for 
any couple of configuration C and C the following rela- 
tion has to hold 



S[C'\C]P 00 (C)=S[C\C'}P 00 (C) 



(46) 



where is the probability distribution we want to sam- 
ple through the Markov chain. 

A simple algorithm to generate configurations so to 
satisfy detailed balance was introduced by Metropolis 36 . 
The idea is to start from a given intial configuration C and 
to propose to visit a new configuration C' with a certain 
transition probability T (C'\ C). Then this new configura- 
tion is accepted or rejected according to the probability 
A (C'\ C) so that the full conditional probability to move 
toward C starting from C is given by 



S[C'\C] =T(C'\C) A(C'\C) 



(47) 



The Metropolis choice for the acceptance probability 
A(C'\C) reads 



A(C'\C) = min\ 1, 



PooCC) T(C\C) 
PooiC) T(C'\C) 



(48) 



It is easy to show that such a choice satisfies the detailed 
balance condition in Eq. (|4"6"1) . While this algorithm is 
completely standard and model independent, two main 
issues have to be taken into account with reference to 
the problem at hand, since they can strongly affect the 
performance or even the reliability of the Monte Carlo 
simulations. 

The first one concerns the choice of the transition prob- 
ability T(C'\C). In the case of interest, we implement 
two main classes of local moves, in which the number of 
kinks in a given channel a is changed by unity, Ak a = ±1. 
These moves amount to add or remove one creation and 
one annihilation fermionic operator in the a channel at 
randomly chosen times along the contour and are re- 
quired for the ergodicity of the matrix S. Indeed it is 



evident that, using these two basic updates any configu- 
ration can be reached, in principle, after a finite number 
of steps 5 ^. Therefore using these two classes of moves 
and the Metropolis acceptance Eq. (J3HJ) we can guarantee 
that the sampling process visits configurations according 
to probability (C). A second issue, which is different 
from the ergodicity one, concerns the efficiency and the 
speed-up of the Monte Carlo sampling (for example the 
number of steps which one has to wait before reaching 
the desired probability distribution). For this purpose, 
additional kind of updates may enhance the sampling 
procedure. In the present algorithm, following standard 
practice in diagMC, we also implement moves that con- 
nect configurations at fixed number of kinks (Ak a = 0) 
such as for example shifting an annihilation operator. 
We also note that other kind of moves, in which more 
than two operators are added/removed/shifted, may be- 
come relevant when dealing with off-diagonal baths. Sim- 
ilarly global moves, in which a whole set of operators is 
changed, has proven to be fundamental in the case of 
multiorbital models 3 -^. We note that, as it happens for 
the imaginary time algorithm^, a major issue in the im- 
plementation of these Monte Carlo moves is to properly 
take into account the structure of the impurity hilbert 
space to avoid, when it is possible, moves toward config- 
urations which have zero weight. In the case of impurity 
models without exchange or hopping terms, these zero- 
weight configurations can be immediately read out since, 
for each channel a, creation and annihilation operators 
have to occur in alternated order along the contour, to 
have a finite local trace in Eq. (|57|) . This leads to a very 
convenient segment picture^. 

A further point which requires some comment concerns 
the evaluation of the acceptance ratio in Eq. (|48| . As 
can be seen from the expression of the weight W (C) 
in Eq. (|57|) . this amounts to evaluate the ratio of two 
determinats as well as the ratio between two local traces. 
While for the former fast update routines are available, 
which makes this operation rather efficient, dealing with 
the ratio among local traces is the most time-consuming 
part of the algorithm, at least in the general case of a 
multi-orbital quantum impurity model. In such a case, 
indeed, one have to evaluate from scratch the whole chain 
of fermionic operators. Several tricks have been proposed 
to implement this evaluation^! in an efficient way and we 
have used most of them in our algorithm. In particular, 
we write the fermionic operators entering in Eq. in 
the basis of local eigenstates and store the whole chain 
of matrix products from left to right (and viceversa), so 
that the evaluation of trial moves is reduced to few matrix 
multiplications. 

In the next section we describe the first application of 
the diagMC algorithm on the Kadanoff-Baym-Keldysh 
contour to the single impurity Anderson Model. In par- 
ticular we will focus on the impurity real-time dynamics 
after a local quantum quench. 



V. REAL-TIME DYNAMICS IN THE 
ANDERSON IMPURITY MODEL AFTER A 
LOCAL QUANTUM QUENCH 

Quantum quenches in strongly correlated systems have 
recently attracted lot of scientific interest, especially in- 
spired by exciting experiments on cold atomic gases 5 
where sudden changes of Hamiltonian parameters has 
been realized and the non equilibrium dynamics moni- 
tored in real-time. In the context of impurity models, in- 
stead, the study of quantum quenches has a long history 
which goes back to the fundamental work by Nozieres and 
De Dominicis on the X-ray edge singularity 3 ^, passing 
through the famous Anderson and Yuval approach to the 
Kondo models. More recently, this problem stimulated 
new interes t 13 ! 40 ' 41 , due to the experimental progresses 
in nanotechnology, which made it possible to contact mi- 
croscopic quantum objects with metallic electrodes, thus 
realizing quantum impurity models in a fully tunable set- 
up 4 ^. Two kinds of quenches can be considered in this 
context, depending on the amount of energy that is in- 
jected into the system, also referred as the work done 
during the quench. Global quantum quenches are partic- 
ularly relevant for transport through correlated nanos- 
tructures, where a net current flow is forced by suddenly 
switching on e.g. a dc bias voltage. Since the switched 
perturbation is extensive, the system is driven into a non- 
equilibrium steady state at long times™. Conversely, lo- 
cal quantum quenches amount to suddenly change the 
impurity Hamiltonian. These kinds of quenches could be 
realized in an optical absorption experiment, as suggested 
in 4 ^, and the resulting non-equilibrium dynamics can 
be tracked in real-time using pump-probe techniques or, 
in real-frequencies, measuring the absorption lineshape. 
Furthermore, local quenches are interesting as they are 
the simplest examples of non-equilibrium processes whose 
statistics may show non trivial fluctuations^ 5 -. 

To test the algorithm, in this section we study the non- 
equilibrium dynamics in the Anderson Impurity Model 46 
after a local quantum quench. This model, which serves 
as a fundamental paradigm for strong correlation physics, 
describes a single interacting fermionic orbital coupled to 
a an equilibrium bath of free conduction electrons. The 
local hamiltonian before and after the quench reads 

n tc [4, c a ) = -y- (n - l) 2 + e d ± n n = ^ c\c a . 

(49) 

The conduction electrons are assumed to be non inter- 
acting, hence the coupling to the impurity occurs via an 
energy-dependent hybridization function r(e), which is 
defined in terms of the conduction density of states (DoS) 
p(e) as 

r(e) = tt l^kl 2 S(e — £k) — ttV 2 p (e) , (50) 

k 

where we have assumed for simplicity Vk independent of 
momentum. As a model for the DoS we start considering 



9 



a flat band of width 2W 



p(e)= Po 6(W~\e\) 



(51) 



which encodes the main physics of a metallic conduction 
bath, namely a finite weight at the Fermi level (po = 
1/2W at e = 0) and a finite bandwidth. In this case 
the basic energy scale describing the coupling between 
the impurity and the bath is the hybridization strength 
r = ttV 2 pq. In all calculations we take T as our unit of 
energy and choose W ~ 1QT. 

This section is structured as follows. We first discuss 
some aspect of the algorithm (statistics of kinks and av- 
erage sign) in the specific case of the Anderson Impurity 
Model and then present the results for its charge and spin 
real-time dynamics after a local quantum quench. 
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A. Performance Analysis of the 
Kadanoff-Baym-Keldysh diagMC algorithm 

In order to analize the performances of the diagMC al- 
gorithm on the Kadanoff-Baym-Keldysh contour we will 
consider two main quantities, namely the probability dis- 
tribution of perturbative orders in the diagrammatic ex- 
pansion and the average sign of the Monte Carlo weight, 
both being sensitive measures to establish the efficiency 
of the method. In the specific case of the single impu- 
rity Anderson Model, with reference to the notation in- 
troduced in section IIII1 we have only J\f — 2 channels, 
corresponding to spin a =t, 4- 



Statistics of Kinks 



As we have shown in section ITVl diagMC amounts to 
stochastically sample the expansion for 2 (t) in powers 
of the hybridization, by performing a random walk in 
the space of diagrams. It is therefore quite natural to 
monitor during the simulation the statistics of different 
perturbative orders (number of kinks), namely the prob- 
ability to visit a Monte Carlo configuration C featuring k 
creation vertices (and k annihilation vertices) in the spin 
sector a. We define this quantity as 



Pa (k) = 



E c \W(C)\5{k{C a )-k) 

E c \w[C)\ 



(52) 



where the Monte Carlo weight W (C) has been defined in 
Eq. (|37ll . The behaviour of this probability distribution 
is plotted, as an example, in the left panel of figure [3] 
for increasing values of the measuring time t* starting 
from = 0, which corresponds to the equilibrium initial 
preparation. We note that all hystograms are strongly 
peaked around an average value k, larger perturbative 
orders k > k appearing with an exponentially small prob- 
ability. Notice that whenever this would not be the case, 
namely if arbitrarily large perturbative orders (k 3> k) 
would give a finite contribution to the result, a diagMC 



Figure 3: Probability distribution of different perturbative 
orders k sampled during the simulation. Data refers to a 
local quench in the Anderson Impurity Model, starting from 
U- = to U+ = 10r at particle-hole symmetry, for T = 0.1F 
and W = 10F. Left panel shows the statistics ok kinks along 
the whole contour, while the top and bottom right panels 
display the hystograms for the Matsubara and Keldysh sector, 
respectively. 



algorithm to work would need an explicit cut-off k max on 
the perturbative order and the result would then require 
an extrapolation 59 to k max —> oo. However for quan- 
tum impurity models, at least for the weak and strong 
coupling algorithm a 21 ' 22 ! 31 , this is not the case. Figure [3J 
confirms that all orders are included and contribute, with 
their own weight, to the final result. This fact ensures 
that the outcome of our diagMC calculation is an unbi- 
ased result which do not correspond to any truncation 
at finite-order of the perturbative expansion but rather 
represents a numerical resummation of a formal expan- 
sion. From figure [3] we note that upon increasing t* the 
whole hystogram shifts toward larger values of k since 
kinks start to be added on the two Keldysh branhces. It 
is therefore interesting to resolve this increased perturba- 
tive order in the two sectors of the simulation, namely the 
Matsubara and the Keldysh one. To this extent, we plot 
in the right panel of the same figure the probability distri- 
bution of having k to t kinks with spin a on the Matsubara 
axis, Pmo- (k) (top right panel), or ktot kinks with spin a 
on the Keldysh branches, Pk<j (k) (bottom right panel), 
where fc tot means that we are considering both creation 
and annihilation vertex. At t* = 0, the initial state, the 
Keldysh branches are empty while the Matsubara sector 
is filled with an even number of kinks (to ensure particle 
conservation). Upon increasing t* > 0, the system starts 
evolving in real-time and we note a transfer of weight 
in the Keldysh sector toward finite values of k. At the 
same time the Matsubara probability distribution does 
not change its center of gravity and rapidly converges to 
a final distribution, now allowing also for an odd num- 
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Figure 4: Scaling of the average number of kinks with maxi- 
mum time t±. Left panel shows data at fixed U+ = 0, tuning 
the strength of the Coulomb repulsion U- = 0, 5, 10. See the 
perfect linear scaling with the same slope a = dk/dt. Right 
panel displays data at fixed U- — 10, U+ = tuning the 
strenght of the conduction bandwidth W . We see how the 
slope a increases with bandwidth making increasingly diffi- 
cult to access large time scales in the regime W S> T. 



Figure 5: Average Sign as a function of time t for different 
initial preparations. We clearly see an exponential decay on 
a vey short-time scale. Left panel shows data obtained fixing 
the final value of the interaction U+ = and tuning the initial 
value U=0, 5, 10. We see a slight increase of the average sign. 
Right panel shows the dependence of f\ from the bandwidth 
of conduction electrons and suggest that much longer time 
scales can be reached in the regime W ~ V. 



ber of kinks (the total particle conservation is ensured by 
kinks in the Keldysh sector) . The scaling of the average 
number of kinks k with measuring time t*, temperature 
and other physical parameters is also relevant and in- 
structive. In the equilibrium case (corresponding here to 
t* = 0), k has been shown in Ref. 31 to be proportional 
to inverse temperature ft with a prefactor given by the 
average hybridization energy per spin, 

h = -miyb)- ( 53 ) 

Since (rihyb) decreases upon increasing the correlation 
strength U the diagMC method in imaginary time works 
extremely well in the regime U ^> T being able to reach 
very low temperatures compared to the energy scales in 
the problem. Unfortunately, the very convenient scal- 
ing of Eq. (|53p does not hold anymore for the real-time 
dynamics, as was also noted in previous works 25 . In 
figure 2] (left panel) we plot k as a function of time t* 
for different initial preparations U- = 0, 5, 10. We note 
an almost perfect linear scaling with time, as expected, 
while the effect of starting from a correlated initial state 
U- 7^ generally helps since it decrease the value of k at 
t = 0. Nonetheless, a finite Coulomb interaction in the 
final state, U+ ^ 0, has no effect on the average num- 
ber of kinks sampled, as shown by the dashed line in left 
panel which exactly lies on top of the U+ = results. 

Summarizing, we conclude that the scaling of the aver- 
age number of diagrams for the real-time algorithm gen- 
erally reads 

K = at + k?, (54) 



a being a costant independent on U. It is therefore 
natural to ask what is the energy scale controlling this 
prefactor. As we show in the right panel of Figure [U 
a strongly increases with the conduction bandwith W 
(and presumably also on the hybridization strength). As 
a consequence of Eq. (]54f) , accessing large time scales in 
the regime W 3> T becomes increasingly difficult with 
this approach. This is due to the fact that both the com- 
putational cost of the algorithm and, in particular, the 
average sign of the MC weights strongly depend on the 
average number of kinks fc, exponentially the former and 
power- law the latter. 



2. Average Sign 

Another important quantity to monitor during the sim- 
ulation is the average sign of the MonteCarlo configura- 
tions, which is tightly related to the accuracy we can get 
on physical quantities at fixed CPU time. Indeed a van- 
ishing average sign turns into very large error bars on 
MonteCarlo averages that makes the simulation unsta- 
ble. In the specific case of the hybridization expansion 
diagMC, it is known that, for what concerns the equi- 
librium (imaginary-time) algorithm, the single impurity 
Anderson Model has always positive signs and that even 
multi-orbital impurity models with rotationally-invariant 
interaction can be efficiently simulated up to moderate 
low temperatures. This situation drastically change when 
dealing with real-time dynamics since even the simple 
non-interacting resonant level model faces a severe sign 
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Figure 6: Quench dynamics in a Resonant Level Model after 
a sudden change of the energy level from Ed- = to Ed+ 7^ 0. 
Dashed lines is the exact solution for n (t) as obtained by a 
standard methods. Points are diagMC results obtained at 
T = 0.1F. We also add the dynamics for the trivial case 
Ed+ = Ed- — to show that unitarity is actually preserved by 
diagMC. 



problem at large enough times scales. This seems to be 
due to the intrinsinc nature of unitary quantum evolu- 
tion and clearly appears from the definition we gave of 
the dynamical partition function Z (t) in section IIIII In- 
deed exact cancellations are built in the whole formalism 
to ensure unitarity. 

In the real-time diagMC one should in general talk 
about average phase, since as we mentioned the Monte 
Carlo weights are complex numbers. Nevertheless this 
quantity, which is defined as 



E c \w(c)\ v (c) 
E c \w(c)\ 



(55) 



turns to be directly related to the probability of visiting 
configurations in the Matsubara sector, namely to the 
probability of having no kinks on the real-time branches. 
As a consequence of this result, which comes directly 
from unitarity, we conclude the the average phase is in- 
deed a purely real number even for t > and therefore, 
without further misunderstanding, we refer to it as the 
average Monte Carlo sign. As we show in figure [5J fj (i) 
depends exponentially on the measuring time t*, namely 
on the length of the Keldysh contour. In the left panel, we 
plot the average sign for different values of the Coulomb 
repulsion C7_ in the initial density matrix. We see that 
fj (t) uniformly increases with U— , namely starting from 
a correlated initial state may result into a larger average 
sign even if the effect is rather small. In the left panel 
we study how the sign depends on the bandwidth W of 
conduction electrons. We see that moving from W = T 
to W = 15r there is a sizeable increasing of the average 
sign, which means that larger time scales can be reached 



Figure 7: Non Equilibrium Dynamics of double occupancy 
D (t) in the Anderson Impurity Model after a local quantum 
quench of the interaction strength at T — 0.1F and particle- 
hole symmetry. In the upper panel we start from an initial 
state with U- = 10F and a very low double occupation and 
quench to different values of U+ /Y = 0, 2.5, 5 from top to bot- 
tom. In the lower panel the opposite protocol is considered, 
namely we start from different values of U-/T — 2, 4, 6 from 
top to bottom and quench to the same final U+=0. In both 
cases we see that after a rather short transient the system 
relaxes to a new equilibrium state. 



with present algorithm in the regime W ~ T. While this 
is not of direct relevance for quantum impurities, it can 
be very intersting for Non Equilibrium DMFT, where the 
conduction bandwidth is of the same order as the cou- 
pling between the impurity and the bath itself. 



B. Charge and Spin Dynamics in the Anderson 
Model after a local quantum quench 

We start by considering the non interacting case, the 
so called resonant level model with [/_ = U+ = 0, which 
allows for an exact solution and can be therefore used 
to benchmark the algorithm. We consider for this sim- 
ple resonant level model a quench of the energy level £d 
that we tune from the on-resonance value e_ = to the 
off-resonance one e + ^ 0. We note that this kind of 
quench can be realized in optical absorption experiments 
with quantum dots, as recently proposed in Ref. H . In 
Figure [6] we plot the real-time dynamics of the impu- 
rity density n (t) for two different quenches, respectively 
above and below the on-resonance value 54 = 0, and com- 
pare the result of diagMC (datapoints) with the exact dy- 
namics which can be obtained using standard methods 47 . 
The excellent agreement with exact results confirms the 
reliability of our numerical approach. 

We then move to the interacting case, namely consider 
a local quantum quench in the Anderson Model with lo- 
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Figure 8: Non equilibrium dynamics of the impurity density 
after a quench of the energy level from Ed- = to Ed+ = 2.5T. 
We compare the dynamics for different values of the Coulomb 
repulsion U — 0, 2.5, 5, 7.5, 10 from bottom to top, revealing 
a much faster thermalization in the correlated case due to 
Coulomb blockade. Inset: thermal value of the impurity den- 
sity as a function of the level position for U = 0, 2.5, 5, 7.5, 10 
from top to bottom. For U S> F the curve is almost fiat 
around Sd — hence departure from equilibrium is suppressed 
after a short transient. 



Figure 9: Spin Dynamics in the Anderson Model after a sud- 
den quench of a local magnetic field. We start from a partially 
polarized impurity, with /i_ = 5.0, T = O.ir and different val- 
ues of local interaction U (top panel). At time t > the mag- 
netic field is switched off and the magnetization is allowed to 
relax toward an unpolarized steady state. We see that upon 
increasing the interaction U the dynamics slows down. Due 
to fine-time resolution we cannot follow the decay of the spin 
toward zero magnetization. However a large magnetic field in 
the final state gives rise to a rather fast relaxation. 



cal Hamiltonian (|4"5)) . In Figure [7] we show the dynam- 
ics of the double occupation D (t) — (rif (t) n± (t)) af- 
ter a sudden quench of the local interaction strength U. 
Two different cases are considered. In the upper panel 
of Figure we start from the same initial preparation, 
U— = 10r, and quench to different final values of the in- 
teraction U+/T = 0,2.5,5 (from top to bottom). In the 
lower panel of the same figure, we start from different 
initial preparations U-/T = 2,4,6 (from top to bottom) 
and quench to the same final state U+ = 0. 

The dynamics at short times, soon after the quench, 
is controlled by the initial density matrix as expected on 
general grounds. After a short time scale, t s hort ~ 0.1/T, 
the system starts feeling the quench and in fact the curves 
in the upper/lower panel start to deviate from/approach 
to each other. The time scale controlling the approach 
to the steady state is set mainly by 1/r - without cou- 
pling to the bath no dynamics for the charge would arise 
at all. However the final value of the interaction also 
affects the dynamics, as one can see from data in the 
top panel of figure [7J We also compare these findings 
with the non-interacting case, where the quench is per- 
formed on the energy level which is suddenly placed out 
of resonance (see lower panel black curve) . In this situa- 
tion the dynamics appears much slower than the previous 
cases, at least a factor of two. In Figure [8] the prob- 
lem of quenching the impurity energy level is considered 
for different values of the Coulomb repulsion U , starting 
from a level which is initially half-filled. As compared to 



the non interacting U = case, the effect of interaction 
is to make the whole relaxation dynamics much faster 
and the steady state value closer to the starting one, re- 
sulting in some sense into a less pronounced deviation 
from equilibrium. This can be rationalized by consid- 
ering how the density depends, in equilibrium, on the 
energy level (see inset): upon increasing the interaction 
the curve n (e c i) becomes flat around Sd = 0, a signature 
of Coulomb blockade phenomenon. As a consequence, 
any perturbation which moves the impurity occupation 
out of integer filling is quickly suppressed on a short-time 
scale. The overall picture confirms what also found in a 
similar investigation with the time-dependent NRG mli~3l. 
namely that charge dynamics is sensitive to high-energy 
scales, thus resulting into a generally fast relaxation. As 
opposite to the charge sector, the dynamics of spin degree 
of freedom is sensitive to the low-energy physics of the 
model. To probe this dynamics in real-time we imagine to 
add a magnetic field h_ to the local Hamiltonian, which 
partially polarizes the impurity, and suddenly switch it 
off for t > 0. The local Hamiltonian Eq. (|4T)j) now reads 

cr 

In Figure [9] we plot the spin dynamics (S z (i)) starting 
from h- = 5r and switching it off, h+ = 0, for differ- 
ent values of the interaction U. Since the final state in 
the absence of Zeeman splitting is fully simmetric, we ex- 
pect to recover, for large enough times, a relaxation to an 
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unpolarized state with (S z ) = 0. We see that this relax- 
ation is very slow and controlled by a time scale which 
increases with U (see top panel), as opposite to what 
found in previous cases, when the dynamics of charge 
degrees of freedom was probed by quenching the inter- 
action or the level position. Indeed, the spin dynamics 
in the strong coupling regime is controlled by the lowest 
energy scale in the problem, namely the Kondo temper- 
ature, as explicitly shown Accessing such a long 
time scale seems so far unfeasable within the present ap- 
proach, since diagMC simulations become increasingly 
unaccurate at large times due to sign problem, as we will 
discuss in the next section. As shown in the bottom panel 
of Figure [9j the effect of a large magnetic field in the fi- 
nal state is to destroy Kondo effect, thus resulting again 
into a fast relaxation toward a new steady state. An 
interesting direction for future investigation is to study 
the spin dynamics in the Kondo regime under the effect 
of more general non equilibrium processes, other than a 
quantum quenches. In this respect the present approach 
can deal with explicitly time dependent phenomena, such 
as for example an oscillating magnetic field, without any 
truncation of the dynamics. 

C. Non Equilibrium Dynamics for a quantum 
impurity in a gapped or pseudogapped fermionic 
reservoir 

In all cases considered above, we observe at large 
times a convergence to a new equilibrium state which is 
the thermal one described by the final hamiltonian 
namely 

(Q(*->°°)) = T ^ [e _^ +] J (57) 

This is explicitly shown by the dashed line in figure [3 
which represent the result of an equilibrium calculation 
done with imaginary time diagMC with the final Hamil- 
tonian H+. We also note that no effective heating arises, 
namely the temperature entering in Eq. (|57p is the same 
as in the initial condition @ . This is due to the fact that 
within diagMC the fermionic reservoir is treated as an in- 
finite system. The onset of thermalization in a quantum 
impurity model is not surprising—, and it is related to 
the fact that the conduction electrons play the role of 
a thermal bath 43 , able to absorbe the energy pumped 
locally after quench, which is dissolved in the interior 
of the bulk. It is worth noting that this feature is not 
generic of any bath - meant as a macroscopic (infinite) 
system - but rather depends on its spectral properties. 
In the present case, as we are going to see, thermaliza- 
tion is related to the gapless nature of the metallic state, 
whose energy spectrum goes down to arbitrarily small en- 
ergies. To further investigate this issue, we consider now 
the out-of-equilibrium dynamics of an Anderson impurity 
coupled to gapped and pseudo-gapped fermionic bath. 
Even though QI models traditionally deal with genuine 



metallic hosts, the problem of gapped (or pseudo-gapped) 
fermionic reservoirs has a vast literature 4 ^—, that re- 
ceived a large boost in recent years. Eminent examples 
of such a physical situation are provided by adatoms in 
graphene sheet or by nanostructures built up with super- 
conducting materials. 

The equilibrium phase diagram of an Anderson impu- 
rity embedded in a non-metallic host it is by now rather 
well establishe d 12 ' 50 ' 53 . As opposite, the non equilibrium 
real-time dynamics in this class of quantum impurity 
models is much less explored. A detailed study of this 
issue goes well beyond the scope of this paper and it will 
be left for future investigations. Here we limit to eluci- 
date the role played by the presence/absence of low en- 
ergy bath spectral weight on the single impurity dynam- 
ics after a local quantum quench. It is worth to notice 
that this issue can be also relevant to study, within Non 
Equilibrium DMFT, the relaxation dynamics of interact- 
ing electrons after quantum quenches. Indeed DMFT 
amounts to solve a quantum impurity self-consistently, 
using the contour ordered impurity Green's function as 
a seed to generate the new fermionic out-of-equilibrum 
bath. 



1. Gapped Fermionic Reservoir 

We start our discussion considering the case of a true 
gapped fermionic bath. In other words we consider as a 
model for the conduction electrons DoS the following 

p (e) { ° ° < |£| < E » (58) 
M ° j -\A) E g <\e\<E g + W ' ^ 

where 2E g is the band gap at the Fermi Level. This den- 
sity of states results into an energy dependent hybridiza- 
tion function T (e) that we define as in Eq. (1501) , namely 
r (e) = nV 2 p(e). A plot of this function for different 
values of E g is given in figure [TD] We note that in the 
following we will adopt as unit of energy the hybridiza- 
tion width r = ttV 2 /2W, the same as in the metallic 
case. 

The equilibrium properties of an Anderson Impurity 
coupled to a gapped reservoir have been studied with 
NRG in—, and more recently with a perturbative ap- 
proach in 5 ^. The model at particle-hole (PH) symmetry 
flows at low temperature to the local moment fixed point 
where the impurity is asymptotically decoupled from the 
bath. Out of PH the model displays a transition between 
local moment fixed point and strong coupling fixed point 
depending on wheter the gap E g is larger or smaller than 
the Kondo temperature. Here we consider for simplic- 
ity the PH symmetric point which correspond to setting 
e c i- = £d+ = in the local hamiltonian Eq. and 
discuss the real-time dynamics for the double occupancy 
D(t) — (n^ (t) 7i|, (t)) on the impurity site after a sudden 
change of the local Coulomb interaction. 

In figure [TU] we show the double occupancy dynamics 
after a quench from U- = 10T to U+ = for different 
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Figure 10: Non equilibrium dynamics for an Anderson Im- 
purity coupled to a gapped fermionic reservoir. We plot the 
real-time dynamics for double occupancy at the impurity site 
after a quench of the interaction from U- = 10V to U+ = 0, 
at particle-hole symmetry Ed+ = £d- ~ and for T = 0.1F. 
Different values of the gap in the bath E g = 0, 1, 2.5, 5.0 are 
considered, resulting into very different dynamics. Contrar- 
ily to the gapless case (E g — 0, black points) which quickly 
approaches the thermal plateau fixed by PH symmetry and 
indicated by an arrow, Dtherm = 1/4, we see that due to the 
finite gap in the spectrum the real-time dynamics slows down 
thus preventing us to conclude on the long time behaviour of 
D(t). However, for very large values of E g (see E g = 5F) the 
dynamics seems actually to reach a steady state where the 
double occupation is different from Dtherm- 



values of E g at T = 0.1F. Due to PH symmetry the 
thermal value of D computed on the final hamiltonian 
H+ has to be equal to D t herm = 1/4 for U+ = 0. In- 
deed, we see that in the metallic case (E g — 0) D(t) 
approaches rather quickly the expected thermal plateau. 
At the same time opening a finite gap E g ^ in the 
bath reflects into a much slower dynamics which pre- 
vents us from a firm conclusion on the asymptotic be- 
haviour of D(t). We note, however, that for large values 
of E g the dynamics seems actually to reach a station- 
ary state which looks quite different from the expected 
one. Such a behaviour could be interpreted in terms of 
the equilibrium properties of the gapped Anderson im- 
purity model which, as we mentioned, at PH symmetry 
flows to the local moment regime with the impurity ef- 
fectively decoupled from the bath. Given such an initial 
condition and taking into account the large value of the 
gap, which strongly affects the bath properties, one can 
rationalize the slowing-down of the impurity dynamics. 
Indeed, in the limit of very large gaps E g — > oo, a free im- 
purity would have no available mechanism to exchange 
energy and relax to the steady state described by H+. 
While this argument could be in principle satisfying to 
explain results plotted in figure [TU1 at least in the large 
gap regime, it does not take into account completely the 



Figure 11: Non equilibrium dynamics for an Anderson Im- 
purity coupled to a gapped fermionic reservoir. We set the 
gap in the spectrum equal to E g — 1.0 and plot the real- 
time dynamics for the double occupancy after a quench of 
the interaction, at T — 0.1T and PH symmetry. Fwo kind 
of processes are considered, namely a quench from U- = 
to U+ 7^ and the reverse process from U- 7^ to U+ = 0, 
which differ among each other for the average work done dur- 
ing the quench, see Eq. (|60|l in the main text. Top panel 
shows data for U- = 10F, U+ = and viceversa, while bot- 
tom panel data for U- = 2.5F, U+ = and viceversa. We 
see that, provided the average work done is sufficiently larger 
than the gap 2E g , a fast thermalization can occurr also in 
a gapped model (see top panel, black points). As opposite, 
when the amount of energy pumped into the system is too 
small the dynamics slows down and we cannot conclude with 
present data wheter thermalization takes place or not. 



nature of a quantum quench process. To further investi- 
gate this point we now reverse the perspective, namely we 
fix the gap Eg in the spectrum and change the strength 
of the quench, namely we change the final value of the 
interaction U+ while keeping fixed {/_ = as well as the 
level position so to ensure PH symmetry. This allows to 
study how the non equilibrium dynamics depends on the 
amount of work done during the quantum quench. As 
was recently suggested ini£ the statistics of the work is 
a key quantity to characterize a non equilibrium process 
such as a quantum quench. Its average value W gives a 
measure of the energy pumped into the system and turns 
to be given—, in the case of an istantaneous quench, by 

W = (H- -H+)-, (59) 

where the average ( • ) _ is taken over the initial equilib- 
rium density matrix p eq oc e~^ u ~ . In the case of a local 
quantum quench such as the one we are considering, the 
average work W is given by 

W = (U--U+){D)^. (60) 

We see therefore that the work W depends not only on 
the strength of the quench, namely the change in the 
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interaction, but also on the initial condition. As we are 
going to see, this quantity greatly affects the resulting 
non equilibrium dynamics. 

In figure [TT] (top panel) we plot the dynamics of double 
occupancy D(t) at PH symmetry after a quench of the lo- 
cal interaction. We set T = 0.1T and choose a fixed value 
of the gap Eg = T. We compare two kind of processes: 
one starting from U- — to U+ ^ and the reverse one, 
which starts from U- ^ and quenches to U+ = 0. Quite 
interestingly we see that, provided the average work W 
is above the threshold of the (semi)gap E g , as for the 
process [/_ = — » U + = 10T (black curve in top panel) 
for which \W\ = 2.5T, a rather fast thermalization can 
occurr also in the gapped model. Notice indeed that the 
expected thermal value for D, which is set by the dashed 
line in figure Qj] (top panel) and which corresponds to the 
value of double occupation computed at equilibrium for 
U+ = 10r, is approached on a rather short time scale. 
We compare these findings with the inverse quench pro- 
cess, starting from U- = IQT and quenching to U+ = 0, 
which in force of Eq. is characterized by a rather 

small average work W < E g . As we see in figure [Til the 
dynamics looks much slower in this case and we cannot 
conclude, on the basis of our data, wheter the thermal 
plateau at D 'therm = 1/4 is actually approached or not at 
longer times. A similar comparative study is performed 
for quenches from U- = to U+ = 2.51" and viceversa 
and the results are plotted in bottom panel of figure [TT1 




Figure 12: Non equilibrium dynamics for an Anderson Im- 
purity coupled to a pseudo-gapped fermionic reservoir. We 
consider the model at PH symmetry and T = 0.1P We fix 
the quench parameters, namely the initial and final value of 
the interaction, equal respectively to U- = 10 and U+ = 0, 
and tune the pseudo-gap exponent from r = (gapless metal- 
lic state) to r = 4. The depletion of low energy states in the 
DoS reflects into a much slower dynamics which eventually, for 
large enough r, seems to prevent the system from reaching the 
value of Dtherm = 1/4 which is fixed by PH symmetry and by 
the choice of U+ = 0. However, due to finite time resolution 
we cannot conclude with present data wheter thermalization 
occurs or not on a longer time scale. 



2. Pseudo-gapped Fermionic Reservoir 



We now consider the dynamics of an Anderson impu- 
rity coupled to a pseudogap reservoir. We consider as 
DoS a pure power-law function, namely 

, s / a\e\ r < |e| < W , . 

^( £ ) = { \e\>W ' (61) 

where a = (r+1) / (2W r+1 ) ensures the proper normaliza- 
tion. This gives rise to a power-law hybridization func- 
tion r(e), that we define in complete analogy with the 
previous cases see Eq. (l50l) . Again, we choose as unit of 
energy the hybridization width T = it V 2 /2W. 

The equilibrium phase diagram of the pseudo-gap An- 
derson Impurity model is extremely rich, featuring at 
particle-hole (PH) symmetry and for < r < 1/2, a 
quantum phase transition at a critical value of the hy- 
bridization T c between a strong coupling regime (for 
r > T c ) where Kondo screening occurs and a local mo- 
ment one (for T < T c ) where the impurity becomes 
asymptotically free at low temperature. As opposite, at 
PH symmetry and for r > 1/2 the only stable fixed point 
is the local moment one and no Kondo effect can be sta- 
bilized for an Anderson Impurity in a gapless reservoir—. 
In the following we will focus for simplicity on this latter 
case (r > 1/2 at PH symmetry) so to avoid any compli- 
cation related to the dynamics across criticality. We note 



that the topic of non equilibrium dynamics across quan- 
tum phase transitions is indeed an extremely intriguing 
problem^, which may deserves furhter investigations in 
the future. However it goes well beyond the scope of the 
present paper. 

As we did for the gapped case, we consider as a starting 
point a quantum quench of the local Coulomb interaction 
between L/_ = 10T and U+ = 0, at T = 0.1T and for 
Ed+ = £d- = 0. In figure [T^] we plot the dynamics of 
double occupancy as a function of time, while tuning the 
exponent r from the metallic case r = to the almost 
gapped one r = 4. As we can see the results we found 
look very similar to what already discussed in the case 
of a gapped reservoir. While for r = 0, namely for a 
finite hybridization at the Fermi level, the dynamics is 
pretty fast in approaching the thermal plateau (we are 
at PH simmetry and U+ = 0, therefore again D t herm = 
1/4), for r ^ the dynamics slows down and for r = 4 
seems actually to get stucked into a non thermal steady 
state. However from these data it is difficult to conclude 
wheter this is really the case or rather that thermalization 
emerges on a very long time scale. 

We conclude this section by discussing how the dy- 
namics in the pseudo-gapped case changes as a function 
of the work done during the quantum quench. To this 
extent we plot in figure Q2] the double occupancy as a 
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Figure 13: Non equilibrium dynamics for an Anderson Im- 
purity coupled to a pseudo-gapped fermionic reservoir. We 
consider the PH simmetric point, at T — 0.1F and for r = 1. 
We study how the dynamics changes while tuning the aver- 
age work W done during the quench (see text). We fix the 
initial interaction U- = and perform a quench to U+ = 10T 
(top panel). We compare the results with the reverse process, 
from U- 7^ to U+ = 0. Despite the power-law DoS, the dy- 
namics in the case of a large quantum quench (large average 
work) turns to be quite fast. On the contrary, in the low work 
regime (bottom panel) the dynamics is much slower and we 
cannot see wheter thermalization take place or not at longer 
time scales. 



function of time, D(t), for T = 0.1T and at PH sym- 
metry. As we have previously done, we fix the initial 
value of the interaction to U- = while tuning the fi- 
nal value U+ (see top of panel) so to change the average 
work W. At the same time, we study the dynamics for 
the reverse process where we fix the final value of the 
repulsion to U+ — 0, while changing the initial condition 
U-. As we have already found in the gapped case, the 
dynamics turns out to be very sensitive to the average 
work done, namely to the amount of energy pushed into 
the system. In particular we can see from figure Q2] that 
quenches with sufficiently large work W can result into 
a rather fast dynamics and thermalization at long times. 
This seems to be the case, for example, of quenches from 
U- = to U+ = 10r (black points, top panel) where 
the thermal value of double occupation with U+ = 10T 
is set by the arrow at the bottom of the panel. In other 
cases, where the work done is not that large, the dynam- 
ics turns to be slow and we cannot conclude about the 
long time behaviour. 



3. Discussion 

To summarize, in this section we have discussed the 
non equilibrium quench dynamics of the Anderson Im- 



purity model in a gapped or pseudo-gapped fermionic 
reservoir after a quantum quench of the local Coulomb 
interaction. For the sake of simplicity we have consid- 
ered only the particle-hole simmetric case and we have 
chosen the parameters in such a way to be always in the 
local moment regime in equilibrium for both gapped and 
pseudo-gapped cases, so to avoid further complications 
due to local quantum criticality which may result into 
very low-energy /long-time scales controlling the dynam- 
ics. 

An important issue we have tried to discuss concerns 
the onset of thermalization at long times. While this is 
expected to occurr for quantum quenches in a conven- 
tional metallic reservoir, one may wonder wheter or not 
the lack of available states close to the Fermi energy could 
result into a lack of thermalization at long times. We 
have shown that the real-time dynamics is strongly af- 
fected, even on short-to-intermediate time scales, by the 
modified spectral function of the bath. In particular, the 
opening of a gap or pseudo-gap at the Fermi level results 
into a slower transient dynamics. While it is tempting to 
explain this fact by invoking the equilibrium properies of 
the model and the mentioned flow to the local moment 
fixed point, one have to take into account also the in- 
trinsic out-of-equilibrium nature of the quantum quench 
process. In this perspective we have identified the work 
W done during the quench as a relevant physical quan- 
tity to describe the non equilibrium dynamics after the 
quench. In particular, for the gapped and pseudogapped 
models, we have shown that a rather fast thermalization 
can occur provided the work done W is sufficiently large 
(for example with respect to the (semi)gap E g ). As oppo- 
site for small quantum quenches characterized by a small 
amount of work done W <C E g , the dynamics turns to 
be much slower and we cannot conclude, with present 
data, wheter thermalization occurs or not on a longer 
time scale, thus leaving the question open to further in- 
vestigations. 



VI. CONCLUSIONS 

In this work we have extended the recently developed 
real-time diagMC method, in its hybridization expan- 
sion formulation, to full Kadanoff-Baym-Keldysh con- 
tour. The resulting algorithm represents a completely 
general and numerically exact approach to real-time dy- 
namics in quantum impurity models, which interpolates 
between the standard equilibrium diagMC 22 defined on 
the Matsubara imaginary-time axis and its recent non 
equilibrium extensions 2 ^— that works on the Keldysh 
contour and requires a special choice of the initial con- 
dition for the dynamics. Merging together these two ap- 
proaches, we are able to deal with the most generic setup, 
namely a strongly correlated quantum impurity model 
initially in thermal equilibrium, which is driven out of 
equilibrium by some external time dependent perturba- 
tion. As a consequence, several kind of initial prepara- 
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tions as well as driving protocols can be considered with 
our approach that allow studying a wide class of non 
equilibrium problems. 

More interestingly, we notice that no constraints are 
required on the nature of the fermionic bath, which en- 
ters in our approach as an input, encoded in the contour- 
ordered hybridization function iAc (t, t'). This allows us 
to deal with the intriguing problem of studying the real- 
time dynamics of the quantum impurity coupled to a 
non equilibrium bath and opens the way to applying our 
method to solve Non Equilibrium Dynamical Mean Field 
Theory. In this perspective the effective action formu- 
lation of the algorithm we have presented at the end of 
section IPPll represents the most natural one. The natural 
extension of this research is to study relaxation dynam- 
ics in correlated macroscopic quantum systems using Non 
Equilibrium Dynamical Mean Field Theory. 

As a first application of our algorithm we have stud- 
ied the real-time in an Anderson Impurity Model after 
a local quantum quench. In the case of a metallic reser- 
voir we have discussed time scales controlling charge and 
spin relaxation. While the former is a rather fast process 
mainly controlled by hybridization T, the latter turns to 
be a much slower process associated with the lowest en- 
ergy scale in the problem, namely the Kondo temperature 
Tjc. As we have shown in section IV B I , the charge time 
scale can be reached within the present approach, while 
the decay of a polarized spin cannot, due to sign problem 
which makes calculations at very long times increasingly 
difficult. 

Finally we have addressed the non equilibrium dynam- 
ics of an Anderson Impurity coupled to a gapped or a 
pseudo-gapped reservoir. Even though we restrict our 
attention to the PH symmetry and to power-law expo- 
nents r > 1/2, for which the equilibrium phase diagram 
in both gapped and pseudo-gapped cases only features 
a local moment fixed point, the real-time dynamics for 
charge degrees of freedom turns out to be rather intrigu- 
ing. In particular we distinguish two regimes depending 
on wheter the amount of work done during the quench, 
W, is large or small with respect to typical energy scale 
in the DoS. In the former case we observe a rather fast 
dynamics which may give rise to thermalization, while in 
the latter case a much slower dynamics which prevent us 
from drawing definite conclusions on the long-time be- 
haviour. The investigation of real-time dynamics in this 
class of quantum impurity models represents, in this per- 
spective, a very intriguing and challenging open problem. 
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Appendix A: Contour-ordered Hybridization 
Function 

In this appendix we discuss with some more details 
the contour-ordered Hybridization Function iAc (t, t') 
we have introduced in the text, which is the basic object 
entering the hybridization expansion on the Kadanoff- 
Baym-Keldysh contour C. This function encodes the ef- 
fect of the bath on the impurity degrees of freedom, as it 
clearly appears in the effective action formulation of the 
theory. As we have shown in section IIII1 in the case of 
a quantum impurity coupled to an equilibrium fermionic 
bath the hybridization function iAc (^1,^2) can be writ- 
ten as 

iA c (h,t 2 ) ee V£( T c fk (h) ft (h)) bath , (Al) 

k 

while in general, i.e. for out of equilibrium fermionic 
baths, a parametrization of this function in terms of 
time independent Anderson impurity Hamiltonian is not 
possible. From the previous expression we see that 
iAc (ii) £2) is given by the contour-ordered bath Green's 
function evaluated at the impurity site. The meaning 
of this function is the following. We consider a bath of 
free fermionic excitations in equilibrium at temperature 
T, whose hamiltonian generally reads 

-Hbath^Y, £ k/kaAa (A2) 
k a 

We take as initial density matrix po the statistical one, 

p — & fibath 

Pin = ~ > ( A 3) 

and define the contour-ordered bath Green's function as 

Sk»(i, t') = -i( T c (fka(t)ft a tf))), t,t'eC 

(A4) 

where both time arguments t and t' live on the three 
branch contour C, while the average is taken over the ini- 
tial density matrix. The contour time ordering Tc acts 
as described in the main text, namely ordering opera- 
tors according to their time argument on the contour C. 
Concerning the contour-time evolution of creation and 
annihilation operators it is defined as usual 

ha(t) = e ffl - f /k«e- H " f . (A5) 

We now discuss the possible time orderings arising 
from the choice of t and t' along the contour C, which 
naturally lead to a 3 x 3 matrix structure for the contour- 
ordered hybridizaton function. 

tA o (t,0 -> iA ab (t,t') a, 6 =1,2, 3 (A6) 

1. Matsubara Sector 

When both arguments live on the imaginary time axis, 
namely t — —ir and t' — —it 1 , we recover the standard 
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Matsubara Green's function, a part from the i-factor 



and 



„33 

3k 



3 a (r,r'H-*<7V (/ka(r)/L(r'))>- (A7) 



We note this function is antiperiodic and time- 
traslational invariant, therefore we can set r' — and 
compute it in the interval r £ [0, 0\. We obtain for the 
hybridization function the standard result used also in 
equilibrium diagMC 



A 33 (r) 



de 



F (e) n (e) e 



(A8) 



where n (e) is the Fermi distribution function and we 
had explicitly introduced the energy-dependent the hy- 
bridization r (e) 

r( £ )=7r^ \ Vk \ 2 5{e - e k ) . (A9) 



2. Keldysh Sector 

When t and t' are both on real-time branches we are 
in the Keldysh subspace. The Green's function acquires 
a 2 x 2 matrix form, depending on the branch position 
(a,b = 1,2) of the two time arguments, and consequently 
also the hybridization function in Eq. (IA1[) can be writ- 
ten as a matrix 



A ab (t,t') 



A 11 (f,f) 
A 21 (f,f) 



A 12 (t,f) 
A 22 (f,f) 



(A10) 



In particular, when t is greater/lesser than t' on the con- 
tour we get the A 21 /A 12 component, which reads respec- 
tively 

A 21 (t,t')=J ^T(e) (l-n F (e)) &-' ie (*"'') (All) 



and 



:( t -f) 



(A12) 



A 12 (t,t') = -f *r(e) n F (e) 



On the contrary when the two times are on the same 
branch we obtain the time-ordered A 11 or anti-time or- 
dered A 12 hybridization function, which reduce to the 
off-diagonal ones depending on the time interval, namely 

A u r*tM -/*>*' a 21 M') rA131 

\t<t' A 12 (t,i') [AU > 



A 22 (t,t') 



t>t' A 12 (t,t') 
t<t' A 21 {t,t') 



(A14) 



We note that, by construction, since we are considering a 
time independent quantum impurity model, all the four 
Keldysh components only depend on time differences t — 
t' . This property is peculiar of a bath which is in thermal 
equilibrium and do not hold in general for other kinds of 
non equilibrium driving protocols or in the case we are 
solving Non-Equilibrium DMFT. 

3. Mixed Sector 

Finally we have to consider the case in which the two 
time arguments live in different sectors. These mixed hy- 
bridization functions usually take into account the short- 
time memory effects, namely the transient correlations 
due to the chosen initial density matrix. We can distin- 
guish two cases, depending on wheter t >c t' or vicev- 
ersa. In the former case, namely when t = —it is on the 
Matsubara branch while t' lies on the Keldysh branches 
we have 



A31 (r,t') 
as well as 

A 32 (T,t>) - 



-T(e) (l-n F (e)) e 

7T 



-e t g — i e t' 



(A15) 



de 



r(e) {l-n F (e)) e 



i € t 



(A16) 

We note that in this case, as usual for off-diagonal terms, 
the contour-time ordering is fixed independently from the 
values of the time arguments since the Keldysh branches 
are always lesser on the contour C than the imaginary 
branch. Moreover, we note there is no difference in the 
value of the hybridization function if the real-time is 
placed on the upper or lower branch, namely 



A (T,i') = A32 (r,f) 



(A17) 



As opposite, when t lies on the Keldysh contour while 



a = 



is on the Matsubara branch we obtain the other 



two mixed components 



A (t,r') = - 



— T(s) n F (e) e 

7T 



T' ^ — i € t 



= A 23 (t,r') 
(A18) 
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